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Abstract 

We prove that when a class of partial differential equations, gener¬ 
alized from the cable equation, is defined on tree graphs, and when the 
inputs are restricted to a spatially discrete, well chosen set of points, 
the Green’s function (GF) formalism can be rewritten to scale as 0(n) 
with the number n of input locations, contrary to the previously re¬ 
ported 0(n 2 ) scaling. We show that the linear scaling can be combined 
with an expansion of the remaining kernels as sums of exponentials, to 
allow efficient simulations of equations from the aforementioned class. 

We furthermore validate this simulation paradigm on models of nerve 
cells and explore its relation with more traditional finite difference ap¬ 
proaches. Situations in which a gain in computational performance is 
expected, are discussed. 

Keywords, partial differential equation, tree graph, Green’s function, 
sparse, simulation 


1 Introduction 


Neurons have extensive morphological ramifications, called dendrites, that 
receive and integrate inputs from other neurons, and then transmit the re¬ 
sult of this integration to the soma, or cell body, where an output in the form 
of action potentials is generated. Dendritic integration is considered a hall- 
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using the cable equation — a one dimensional reaction-diffusion equation 
that governs the evolution of the membrane potential V(x,t), and is de¬ 
fined on a tree graph representing the dendritic arborization. Inputs, such 
as synaptic currents, are usually concentrated in a spatially discrete set of 
points and depend on this potential through their driving force (Kuhn et ah, 
2004), and possibly through other non-linear conductances such as NMDA 


(Jahr and Stevens, 1990). 


Traditionally, this integration is modeled using compartmental simula¬ 
tions, where space is discretized in a number of compartments of a certain 
length and a second order finite difference approximation is used to model 


the longitudinal currents (Hines, 1984). As the error of this approximation 


depends on the distance step, complex neural structures require many com¬ 
partments and are computationally costly to simulate. 


Often however, stretches of neural fiber behave approximately linear (Scho 


sn 


et ah, 2012), and one is not interested in the explicit voltage at all locations, 


but only at a specific output location. For this reason we proposed the idea 
of using the Green’s function formalism to simulate neuron models that re¬ 
ceive inputs at a limited number of locations (Wybo et ah, 2013). Given a 


number of n inputs, the voltage at an output location Xi can be written in 
this formalism as: 


n j-t 

V(xi,t) = Yl / g(xi,Xj,t- s)Ij(s)ds. 

3=1 J ° 


( 1 ) 


Problems arise when these inputs depend on the local voltage. Since, in such 
a case, all local voltages have to be known, a system of n Volterra integral 
equations has to be integrated: 


n ft 

V(xi,t) = J2 g{xi,Xj,t - s)I j (s,V(x j ,s))ds, i = 1,... ,n. 

3=1 


( 2 ) 


It can be seen that this system contains n 2 convolutions. This unfavor¬ 
able scaling, along with the fact that the convolutions themselves are costly 
to compute and the restriction to point-source non-linearities, significantly 
impedes the computational efficiency of the GF formalism, and restricts its 


usefulness to very small numbers of input locations (Koch, 1998, page 59-60) 


In this work, we are able to significantly improve computational efficiency 
compared to the classic GF formalism by showing that all three perceived 
disadvantages can be overcome. Using a transitivity property for the Green’s 
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function (Koch, 1998) (see appendix qAl) , we show that, when the input 


locations are well chosen, a transformation of the system ([2]) exists so that 
only 0(n ) kernels are required. We term this the sparse Green’s function 
(SGF) formalism. As an example of how an efficient integration algorithm 
can be designed for the resulting system of Volterra integral equations, we 
show that the kernels can be expressed as sums of exponentials using the 
vector fitting (VF) algorithm (Gustavsen and Semiyen, 1999) (see appendix 


and that consequently the convolutions can be computed recursively. 
Finally, we show (in a simplified setting) that when the spacing between the 
input locations becomes sufficiently small, the SGF formalism reduces to the 
second order finite difference approximation (in the spatial component) of 
the original equation. As a consequence, the SGF formalism can be seen 
as a ‘generalization’ of the second order finite difference approximation to 
arbitrary distance steps, as long as what lies in between the distance steps 
is approximately linear. 

We validate this novel SGF formalism by reproducing two canonical re¬ 


sults in neuroscientific modeling. First, we reproduce the result of (Moore 


et ah, 1978) on axonal action potential velocity. Second, we compare our 


SGF formalism with the de facto standard NEURON simulator (Carnevale 


and Hines, 2006) in the case of dendritic integration with conductance based 


synapses. In a final section, we discuss in which cases the SGF formalism may 
yield computational advantage over canonical second order finite difference 
approaches (of which the NEURON simulator is an example). 


The system of equations 

Each edge of the tree graph represents a segment of the dendritic tree, for 
which the cable equation has the following form: 

dV it a 2 d 2 V n 

2wac m -g-(x,t)+—-g- i (x,t)-2irag m V(x,t)+2_,Ic(x,t) = 2^Ii{t,V{x u t))d{x-Xi), 

(3) 

where Cm,g m and r a denote, respectively, the membrane capacitance, the 
membrane conductance and the axial resistance, a denotes the radius of the 
dendritic branch, I c the current contribution of a channel type c and the 
input current at location X{. The ion channel current can depend non-linearly 
on the voltage and a number of state-variables: 

h{x, t) = f c (V(x, t), y c (x, t)), (4) 
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where the state-variables y c (x,t) evolve according to: 


ycj(x,t) 


y c ,j,inf(,V(x,t)) -y c>j (x,t) 
T c,j(V(x,t )) 


( 5 ) 


with T c j and y c ,j,mt functions that depend on the channel type. Lineariz¬ 
ing these currents yields a quasi-active description (Koch, 1998) of the ion 
channels: 

ic,linOM) = ^V(x,t) + '52-7^ S -y e j(x,t), (6) 


dye 


with 


MM) = w 


d ( yc,j ,inf 


T, 


V(x,t) - y c ,j(x,t ), 


cj 


T, 


(7) 


CJ 


where all derivatives, as well as r c j, are evaluated at the equilibrium values 
of the state variables. If there are a total of K state-variables associated 
with ion channels, a system of K + 1 PDE’s of first degree in the temporal 
coordinate is obtained, which can be recast into a single PDE of degree K+ 1. 

Consequently, we are interested in the GF of PDE’s of the following form: 


L{x)V(x, t ) = 


~ <9 2 d 

Lo ^dx 2+Ll ^dx + L2 ^ 


V(x,t) = L 3 (x)Y^ Ii(t)5(x-Xi), 


i =1 


(8) 

where Li(x),i = 0,..., 3 are operators of the form L^x) = Y.k=o Ckyjk, 
K G N, and S is the Dirac-delta function. We assume: 


(i) that an equation of the form (J8j) is defined on each edge of the tree 
graph (let E denote the set of edges). 

(ii) that on each leaf (let A denote the set of leafs) a boundary condition 
of the following form holds: 


Li-^ vex (t) + L x 2 V^(t) = L x 3 l\t ) VA G A, 


(9) 


where L x are operators defined analogously to Li{x ) and where V ex is 
the field value on the adjacent edge e\ in the limit of x ex approaching 
the leaf. Note that in this general form, equation (J9]) can represent 
sealed end (L x = 1, = L 3 = 0) or voltage clamp (L\ = 0, L 2 = 

L 3 — 1 and I(t) constant) boundary conditions or, when there is only 
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one neurite leaving the soma, a lumped soma boundary condition (see 
(Tuckwell 1988), and where the operators can possibly contain higher 
order derivatives if quasi-active channels are present). 


(in) that at each node that is not a leaf (let <F denote the set of nodes that 
are not leafs, and let E(4>) denote the set of edges that join at node 

4 > e <h): 


E l e £-V e {t) = LiV*(t) + L*I+(t), V0 

emo) dx 


e 


where the operators are again defined as above and V e denotes 
the voltage on edge e in the limit of x e approaching the node </>. The 
first condition then expresses the continuity of the voltage at a node, 
whereas the second condition can signify the conservation of current 
flow ( L e = 7T a 2 e / r e a , Lf = Lf = 0) or a somatic boundary condition 


when multiple neurites join at the soma (see again (Tuckwell, 1988)). 


Algorithms to compute the GF of this system of PDE’s have been described 
extensively in the neuroscientific literature: the algorithm by (Koch and 


Poggio, 1985) computes the Green’s function exactly in the Fourier domain 


whereas the ‘sum-over-trips’ approach pioneered by (Abbott et ah, 1991) 


(see ( 

Bressloff and Coombes, 

1997) for another overview) and extended by 

(Coombes et ah, 2007; 

Caudron et al. 

2012 

) uses a path integral formalism. 

We implemented the algorithm given 

in ( 

Koch and Poggio 

1985 

) as we are 


interested in the GF in the Fourier domain. 


2 Methods 

A sparse reformulation of the Green’s function formal¬ 
ism 

In this section we prove formally that a sparse reformulation of equation ([2]) 
exists. Fourier transforming this equation gives the GF formalism in the 
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Fourier domain (Butz and Cowan 1974 Wybo et al. 2013), : 


1 Fi(w) ^ 

( 

v KM ) 

V 


011 M ••• 9ln(u) ^ 
^ fiVtlK) ' ' ' Qnnijd') j 


< h(co) > 

y KK) J 


(ii) 


where Vi,... ,V n denote the held values resulting from the inputs Ii,... ,I n 
at locations 1 ,... ,n (which may be arbitrarily distributed along the edges of 
the tree graph), and (jij denotes the GF evaluated between points i and j (see 
(Koch and Poggio, 1985) for an example of an algorithm to evaluate these 
functions). Note that for the GF, we dropped the ~ that signifies the Fourier 
transform for notational simplicity. Whenever the argument of a kernel is 
not explicitly mentioned, it will be implied that it is the Fourier transform 
that is under consideration. Note furthermore that the GF between points 
i and j is equal to the GF between points j and i, so that the matrix in 
equation ( JIl| is symmetric (Koch, 1998, page 63). Note finally that while 
the input current A may depend on the local voltage (cf. equation (|3j)) , it 
is still possible to take the Fourier transform by considering I i (t,V{x i ,t)) as 
an a priori unknown function of time (= /,;(£)). 


We can rewrite the set of equations (11) so that the held at one location 


depends on the input only at that location and the held at all other locations: 

Vi{u) = fi(u)ii(u) + £ (12) 

where (with G denoting the matrix for which (£■ = and G -1 its matrix 
inverse and omitting the argument c o for notational clarity) 


fi = 

hi, = —{G~ 1 ) i j/(G~ 1 ) 


(13) 


Intuitively, the held at any location can only depend on the local input 
and the helds at the neighboring locations. After introducing the necessary 
dehnitions and notations, we will prove this intuition formally. 

Notation 1. Let A denote a matrix. A[j\ i] then denotes the same matrix 
with the j ’th row and i ’th column deleted. 
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This way, (—l)- 7 " 1 "* det(H[j; i]) gives the (j, i)’th minor of A and the ele¬ 
ments of the adjugate matrix of A can be written as: 


adj(^4)ij = (-l) J+ *det(j4[j;i]) 


(14) 


Then, by using Cramer’s rule (see for instance (Horn and Johnson 2012, 
pages 17-24)): 


A~ x — 


adj(A) 
det(H) ’ 


(15) 


equation (13) can be written by as: 


fi = det(G)/ det(G[f; *]) 

hij = (—1)* +J+1 det(G[j; i])/det(G[i; i]). (16) 


We specify the discrete set of input locations as follows: 

Notation 2. V denotes a set of n points distributed on a tree graph. 

The geometry of the tree graph will induce a nearest neighbor relation 
on V. We define this relation in the following way: 

Definition 1. Two points i, j G V are nearest neighbors if no other point 
lies on the shortest path between them. 

This allows for the definition of sets of nearest neighbors: 

Definition 2. A set of nearest neighbors J\f is a subset of V in which each 
pair of points is a pair of nearest neighbors, and for which no other point can 
be found in V that is not in J\f but is still a nearest neighbor of every point 
in Af. 


In Fig[l]4 we show an illustration of these sets. 


Notation 3. We use m to denote the number of sets of nearest neighbors 
which can be found in a given V. 


Definition 3. For any set of points V = {1, ... ,n}, G(fP) is the matrix of 
transfer kernels: 


GiV) 


1 fill 

512 ' 

filn 

fifiii 

5Vi2 ' 

Qnn J 


(17) 
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Figure 1: Schematic of the rationale behind the sparsihcation of the Green’s 
function formalism for tree structures. A: Illustration of the sets of nearest 
neighbors J\f q (q — 1,2) induced by the tree structure. For any pair of points 
in A f q there is no other point on the shortest path between them. B: Illustra¬ 
tion of how the structure of the matrices A and A[j, i\ (here j < i) gives rise 
to the sparseness of the SGF. These matrices can be written in ‘block upper 
triangular form’ (white: zeros, black: blocks on the diagonal, grey: other 
non-zero elements). The diagonal of A[j, i] (red) is shifted between the j’th 
row and i’th column (blue) compared to the diagonal of A. When j and i 
do not belong to the same set of nearest neighbors, there is always at least 
one zero on the diagonal of the resulting matrix (here, there are two zeros — 
green circles). 




















































However, when V denotes a set of input locations distributed on the 
tree graph, we will only make the argument of G explicit when we consider 
a subset of V. Hence it is understood that G = G(V), as is the case in 


formulas (13) and (16). 


Definition 4. An attenuation function a t] is defined by: 

a ij = 9ij / 9jj ■ 


(18) 


Note that trivially an — 1. 


Definition 5. AifP) is the matrix of attenuation functions a % j between any 
two points i,j G V (analogous to definition^. 


The significance of these sets J\f (definition [2]) is due to the the following 
lemma (Koch, 1998, page 63): 


Lemma 1. Transitivity property. If i , j E V are not nearest neighbors, 
then for every point l on the direct path between them it holds that: 


9ij 


9u9ij 

9u 


(19) 


Proof. See appendix §|Aj 


□ 


This leads directly to a similar transitivity property for attenuations: 


(/. /j (iij. ( 20 ) 

Note that we may obtain the matrix A (definition [5]) from G by divid¬ 
ing, for all z = 1,... ,n, its Tth column by g vi , and that as a consequence 
(n, 9ii) det(H) = det(Cr). Hence (16) can be written as: 


fi = gn det(H)/det(H[i;i]) 

hij = (—l) l+ - 7+1 det(H[j; i])/ det(H[i; i]). (21) 


Given these definitions and notation, we can capture the aforementioned 
intuition formally as follows: 

Theorem 1. Consider an arbitrary set of points V on a tree graph. Then 
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i) for a point i that is in exactly p sets A/"i,... ,Af p , 

det(A(A/"i)) • • ■ det(A(A/’ p )) 


fi Qi\ 


det(A(A/”i U ... U A fp)[i]i}) 


( 22 ) 


ii) for a pair of points i,j that are nearest neighbors, and for which there 
consequently exists a set N 3 i,j, and where there are exactly p other 
sets Af\,.. ., A f p that contain i: 

h - ( i y+. 7 +i det ( A ( A/ i)) • • .det(A{Af p )) det(A(AQ[j;i]) 
i3 1 J det(A(A/”i U ... U A/" p UA/")[i;i]) 1 J 

Hi) for a pair of points i,j that are not nearest neighbors, 


hij = 0 and hji = 0. 


(24) 


Remark 1. To unclutter the notation, we use the indices i, j G A f to refer 
to the actual points in V, their corresponding positions in the full matrices 
A and in the restricted matrices A(Af). This is justified, as it does not 


influence formulas (22) and (23). Permuting the numbering of a pair of points 


amounts to permuting the corresponding rows and columns of the matrices, 
and thus does not change the determinants. Restricting A to A(Af) amounts 
to deleting the appropriate rows and corresponding columns from A, and 
thus does not change the factor (—l)*+-? +1 . 

Proof. For convenience, we introduce the following ordering scheme for the 
points: we choose one point as the root of the tree and give it the lowest 
index, and then reorder the other points such that within each set A f, the 
point closest to the root has the lowest index, and the other points in A f are 
numbered consecutively. We then reorder the sets A f so that the set with 
highest index contains the point with highest index, and so on. 

Let there be a total of m sets of nearest neighbors A f q , q = l,...,minP. 
First, we show that 

m 

det(A) = ](([ det(A(A/”q)) (25) 

9=1 

by applying elementary row operations recursively set after set. Let’s start 
with A f m , the last set of nearest neighbors, containing k elements. Our num¬ 
bering scheme guarantees that the last k — 1 rows of A correspond to the 
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k — 1 elements in Af m that are not closest to the root. Let us denote the 
point in Af m that is closest to the root by l. Between every d ^ Af m and 
every e E Af m the transitivity property ( |20j ) holds, and thus a e d = a e iaid . 

By applying row operations R e (A ) —» R e (A) — a e iRi(A), that do not 


change det(H) (Horn and Johnson, 2012, pages 9-10), and using (20) along 


with the trivial identity a e i = a e idu , the last k — 1 rows of A become: 


R,(A) = 

( 0 ... 0 CL e n —]^- 1_2 CLqiCLi^—J ^-|_2 • • • &en ^el^ln ^ • ( 26 ) 

Thus, the determinant reduces to 

A g A “) =det(^')det(B), (27) 

where A' and A" are the parts of A unchanged by the row operations, and 
B consists of the non-zero elements of the part of A affected by the row 
operations. It holds that 


det(H) = det 


det(.B) = det(H(.A/" m )). 


(28) 


Construct A(Af m ) by taking the last k— 1 rows and columns from A, and the 
attenuations to (from) point l as the first row (column) of A(Af m ). Then, by 
applying the row-operations R e (A(Af m )) —>■ R e (A(Af m )) — a e iRi(A(Af m )), the 
matrix B will be found as the only non-zero minor of the first row elements 
- the determinant of which is multiplied by 1 to give the determinant of 
A(Af rn ). Thus we have factorized det (A) in det(H') det(H(A/" m )). 

Our numbering scheme for the points guarantees that a similar reduction 
can be applied to det(H'), as long as it contains distinct sets Af, which by 
induction on m proves (25). Note that we could have achieved the same 


reduction by applying the column operations C e (A) -E C e (A) — ai e Ci(A ) in 


an analogous manner. Then the matrix in (27) would have its zero part in 
the upper right corner. 

Next we compute the determinant of By removing both the z’th 

row and column, it is as if i is removed from the set V completely. Conse¬ 
quently, Af\ U ... U Af p \ {?} forms a new set of nearest neighbors, leading to 
a factor det(H(A/"i U ... U Af p \ {*})) = det(H(A/"i U ... U Af p )[i; i]) instead of 
det(H(A/i)) ■ • • det (A (A/),)) in the final product (25). This proves expression 


( 22 ). 
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Now consider a pair of points i , j as in point ii) of the theorem, ff either 
i or j is the point closest in Af to the root, we reorder the points in V as 
described above, but now with the extra constraint that within the set Af, 
all points are numbered consecutively. Note that it is always possible to find 
such an ordering for any one set Af (but not for all sets at the same time). 
Assume furthermore that in the ordering of sets, Af comes at the fc’th place 
(k < m). Again, we start by factorizing out the determinant of A(Af m ). If 
k < m and if /, the point in Af m closest to the root, is not j, the factorization 
can be carried out as explained above by using the row operations R e (A ) —>- 
R e (A ) — a e iRi(A ) for e G Af m . If / = j, the deletion of the corresponding 
row prevents us from using it to execute these row operations, but the j’th 
column can be used instead: C e (A ) —> C e (A ) — ai e Ci(A). Induction on m, 
until the set Af is reached, then factors out the determinants of the matrices 
A(Afq) , q > k. If k = m, the previous induction can be skipped. 

When i nor j are the point closest to the root in set Af, further factor¬ 
ization can proceed unhindered by using row operation, leading to a factoriza¬ 


tion similar to (25), except for the replacement det(A(A4)) —» det(A(A4)[j; *]). 

When i (resp. j ) is the point closest to the root, our special numbering 
scheme allows the application of Rd —^ R<i — ddiRi (resp. Cd —t Cd — aidCA) 
for d G Af\ U ... U Afc-i, resulting in the matrix: 


det(A(A/”i U ... U Afk)) = det 


0 


resp. 


det(A(A/i U ... U Af k )) = det 


A(U)[j-,i] 


A' 1 


= det (B) det(A(W)[j;f]), 

(29) 


A(Af)[j-i] J = det (5) det (A(Ar)[j;f]). 

(30) 

It can be checked that in both cases det (B) = det(A(A/"i U ... U A4~i))- 


Further factorization can then proceed unhindered, leading again to (25) 


with det(A(A/l.)) —> det(A(A4)[j;i]). This proves expression (23). 

Finally, consider a pair of points i > j that are not nearest neighbors. We 
assure that the points are numbered in such a way that all the sets Af that 
contain j have points with indices smaller than i (note that since i and j are 
not nearest neighbors, none of these sets contains i). The familiar reduction 
scheme for A[j\ i] by applying row operations can be thought of as writing 
this matrix in ‘block-upper triangular’ form. Its determinant is then the 
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product of all determinants of the block-matrices on the diagonal. Between 
the j'til row and i’th column of this matrix, the diagonal is shifted by 1 to a 
lower row compared to the diagonal of A. Since % does certainly not belong 
to the same ‘block’ as j, there is at least one 0 on this diagonal (Fig [ljB) . 
This proves that hij = 0. That hji = 0, can be proved analogously, but 
now the matrix is written in ‘block-lower triangular’ form by using column 
operations. As such, the statement (24) is proven. □ 


Corollary 1. Number of kernels. When a total of m sets of nearest 
neighbors J\T \,... ,Af m exists within a set of n points V, the number of non¬ 


zero kernels M in equation (12) is: 


M — n + ^ IA f q ( A f q 


9=1 


- 1 


(31) 


with 


A f q the cardinality of the set A f q . 


Proof. First, we prove equation (31). For n points, there are n kernels /*. 


In every set A f q , there is kernels h t] between for every combination of points 
i, j G A f q , with i j. Consequently, within a set A f q , there are A f q (A f q — 1) 
kernels. Hence, in total, there are n + A f q (A f q — 1) kernels. Q 

Corollary 2. Sparseness. For a given tree graph and a given positive 
number n, the configurations of n points that minimize the number M of 


kernels in (12) have: 


M = 3n — 2. 


(32) 


Proof. When n = 1 the statement is trivial. For an optimal configuration 
with n > 1, A f q = 2 for all m sets Af qi and m = n — 1. Hence M = 
n + 2 (n - 1) = 3 n - 2. □ 


From the viewpoint of computational efficacy, corollaries [I] and [2] indicate 
that there are ‘well-chosen’ and ‘badly-chosen’ ways for the input locations to 
be distributed on the tree graph. For a ‘well-chosen’ configuration |A/"| = 2, 
or at least |A/"| -C n, for all A f, whereas in worst case scenarios there is only 
a single set A f. 
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An efficient method to integrate the system of Volterra 
integral equations 

Transforming ( []~2| ) to the time domain results in a system of Volterra integral 
equations: 

Vi(t)=f fi(t-s)I i (s,V(x i ,s))ds + '52 f hij(t • s)Vj(s)ds Vi, (33) 

JO ■ / ; j 0 

where all kernels h t j between points i,j that are not nearest neighbors are 
zero due to theorem [lj In this form, the SGF-formalism is well-suited to 
simulate neuron models. 

Let O(rifc) denote the number of operations required to compute a convo¬ 
lution each time-step. If these convolutions were to be integrated naively, by 
evaluating the quadrature explicitly (the Quad approach), n k would denote 
the number of time-steps after which the kernel can be truncated. However, 
the kernels are computed in the frequency domain, so if a partial fraction 
decomposition in this domain can be found: 


Li ryl 

m « E - 


T^i iuj + a i 


lij 


Hij{w) ~ ■ , i ) 

t^l luJ + a ij 


(34) 


where 9ft(a) < 0, then in the time domain the kernels can be readily expressed 
as sums of one sided exponentials: 


MO ~ 


l=i 

Lin 


(35) 


M*) K 


1=1 


when t > 0, and fj{t) = hp-(f) = 0 otherwise. Such a decomposition can be 


derived accurately by employing the vector fitting (VF) algorithm (Gustavsen 


and Semiyen 1999) (see appendix qBJ) . Consequently, the convolutions with 


the kernels become sums of convolutions with exponentials, which can be 


readily expressed as simple differential equation (see for instance (Lubich 


and Schadle 2002)). We call this the Exp approach and n* is the number 
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of exponentials per kernel in this case. Usually, for the Exp approach, n& is 
smaller than for the Quad approach. Nevertheless, many of the exponentials 
of the VF algorithm are used to approximate the small t behavior of the 
kernels, and hence have a very short time-scale. The large t behaviour is 
often described by one or a few exponentials. This suggest that a ‘mixed’ 
approach could yield optimal performance, where for small t the quadrature 
is computed explicitly, and for large t the convolution is computed as an 
ODE (or the sum of a few ODE’s). We now give a detailed account of the 
mixed approach. 

Let t = Nh, with h the integration step and N a natural number, and let 
us assume that V^r) is known for all i and for r G {t — kh ; k = 1 ,,N}, 
and that we want to know Vt(t + h), with h > 0. We split the convolutions 
in equation ([33]) into an quadrature term (/XX) an d a exponential term 

(fo~ Kh ) : 

pt-\-h _ pt~\~h 

Vi(t + h) — / fi(t + h - s)Ii(s)ds + X] / + h- s)V 7 (s)ds+ 

Jt—Kh Jt—Kh 

rt—Kh _ rt—Kh 

/ fi(t + h - s)Ii(s)ds + X / hij(t + h - s)Vj(s)ds, 

J 0 . / _• J Q 

(36) 

where K is a natural number that needs to be chosen. We assume that 
between the temporal grid points t — kh and t — (k + 1 )h, Vj(t) can be 
approximated by a linear interpolation: 

Vj(s) =3 V j (t-kh)+ V ’ i ' t ~ ~ V ii l ~ ( k + 1 ) ft ) ( a -t+kh), t-(k+l)h < s < t-kh. 

(37) 

With this and the exponential approximation (35), the quadrature term for 
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convolutions with Vj becomes: 


rt-\-h 

/ hijit + h — s)Vy(s)ds ~ 
Jt-Kh 

W + ft)£ 


f vj (t-kh)Y, 

k =0 l 

+ Vj (t-Kh)Y ; 


,2k + J 

<4jh 


Hj I „ot[.h _ Y 


e ^ 


7 


a\]h 


e a l i:j {k+l)h _ 2 e a ij kh _). 


Aij „a‘..Kh , 7i' 


-p^ij -L 

/ e ^ / 2 l 

«ii a ij h 


3? 


(38)' 

and similarly for the convolutions with the input /*. For the exponential 
term, we get the following: 


rt—Kh , pt—tin , 

/ /^(t + h - s)V j (s)ds = 7 \ je a ^ K+ ^ h / e a y (t-A-fc--)^.( s ) dSj 

0 7 ^ 0 

(39) 

where — Kh) = Jq~ K h 7 ^ e Q o ■ /- ~ (s) ds is the solution of an initial 

value problem: 

uUt) = a^uUt) + lljVjit) 


t-hK 


<i(t = o) = o, 


(40) 


whose value can be computed recursively: 


, rt—Kh , 

u\At - Kh) = e< h uUt - (K + l)h) +/ e a ^ s V j (s)ds. (41) 

J J Jt-(K+l)h 


~{K+l)h 


Using again the linear approximation (37), this becomes: 


uUt - Kh) « e a v h uUt - (K + 1 )h) 


o' 


+Vj(t-Kh) 


7 1 


,% + 

a ij a ijh 


ij ( a 1 . h 

I p ij 


(42) 


+V/(t - (I< + 1 )h) 


Aij a 1 - h Aij ( a l ..h -i 
—r-e «- Kj- e o — 1 

«ij h \ 


Analogous considerations apply for the convolutions with the input /j. 
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Let us now define a matrix H 0 by: 


(#»)« = E 


Jij_ lij 

1 T 


a 


u 


h 


e< h 


(43) 


and a tensor Hi by: 


= E 

i 

W)S = E 


_ 2 e a ij kh + e"u (fc_1)hN ) 

Q!jj h V ) 


lij at kh , Tij / a 1 , .kh a l ..(K—l)h 

— f~e « H- ToV ' e lJ ~ e aO ’ 

qlaH~l 


a 


v 


y 


for fc = 0,..., A' — 1 

for k = K. 

(44) 


when i, j are nearest neighbors, and (Hq)^ = (iTi)L = 0 otherwise. For the 


input kernels /), we define a vector F 0 (with (F 0 )j an analogous to (43)) and 
a matrix F] (with (F^ analogous to (44)). Using these definitions, we may 


write equation ( |36| ) as: 

^i(t + /i) = 




(W + ft) + E(^i)W^^) 


/c=0 


E 


K 


(Ho)<jV,(t + ft) + E( ff i)fU(* - kh) 


k =0 


(45) 


+ y - Aft) + E E e“ , « (K+1) \4,(1 - Aft). 

I j^i l 

To simplify the notation, we group all terms that do not contain voltage or 
input at time t + h in a vector k(t), for which: 


ki(t) = 


K 


E (F^hd-kh) 


k =0 


+ E 


K 


E(Hi)vW ~ W 


k =0 


+ £ e< {K+l)h u\(t - Kh) + ]T ]T e a ^ K+1) %(t - Kh) 

l j^i l 


(46) 


Consequently, equation (45) becomes in matrix form: 


(I - H 0 )V{t + h) = diag(F„)I(t + h) + k (t). 

Solving this matrix equation then gives V(£ + h). 


(47) 
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Note that when the input current Jj depends on the local voltage, equation 
(47) is only semi-implicit, since the voltage at time t would be required to 
compute the current at time t + h. This description may be unstable for 
certain ion channels. However, the direct dependence of most currents on 
the voltage is linear (in the case of ion-channels for instance, non-linearities 
arise through the non-linear dependence of state variables on the voltage): 


Ii(t,V(xi,t )) = c-iit) + di{t)V{xi,t). 


(48) 


Using this, equation (47) can be modified in the following way: 


(I - H 0 - diag(F 0 © d (t + h))V(t + h) = diag(F 0 )c(t + h) + k(t). (49) 

where © denotes the element wise multiplication. This way of integration, 
while still being semi-implicit (since the equations for state variables of ion 
channels (|5]) are still integrated explicitly) is stable and standard in neuro¬ 
scientific applications (Hines, 1984). 

Note furthermore that the matrix (I — Hq) resp. I — i/o — diag(Fo©d(f + /r) 
is structurally symmetric: whenever an off-diagonal element is nonzero, its 
counterpart opposite to the diagonal is nonzero as well. When \J\f\ = 2 for all 
A/", and when the input locations are ordered in the right way, this matrix is 


a Hines matrix 

Hines 

198' 

1). For such a matrix a linear system of equations 

of the form (47 

resp. 

(49) 

can be solved for V in 0(n) steps instead of the 


usual 0(n 3 ) steps. 

Note also that we did not discuss the initialization steps of this algorithm 
explicitly. We omitted this discussion since, in neuroscience, it usually is 
simply assumed that the neuron model is at equilibrium at all times t < 0. 
Hence, the vectors k(t = 0), k(f = h), k(t = 2 h ),... can be computed easily 
by assuming that V(—kh ) = 0 for all required k's. 

Finally, we remark that K, the parameter in the algorithm which deter¬ 
mines the limit Kh below which the quadrature is evaluated explicitly, can 
be chosen to minimize the number rik of operations per kernel. In principle, 
K could be chosen separately for each kernel. In the present derivation we 
however opted not to do so for simplicity. 


Nonlinear terms and the small Ax limit of the SGF for¬ 
malism 

While a set of ion channels that behaves approximately linear may be incor¬ 
porated directly in the SGF formalism, the question what to do with channels 
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that act truly non-linear remains. Such channels can be moved from the left 
hand side of equation ([3]) to the right hand side, and treated as an input 
current that depends on the local voltage. This current needs to be ‘com¬ 
partmentalized’, i.e. expressed at a discrete number of input locations, with 
a certain separation Ax. 

A ‘recompartmentalizatioiT of this type leads to the question of what the 
relation between the SGF formalism and the canonical second order finite 
difference approximation is. In both approaches, the input is compartmen¬ 
talized. In the finite difference approximation however, the entire cable (i.e. 
the longitudinal, capacitive and leak currents) is compartmentalized as well, 
whereas in the SGF formalism, the cable is treated exactly. Consequently, in 
spatial regions with truly non-linear ion channels, the accuracy of the SGF 
is equivalent to (or better than, since the linear currents are still treated 
exactly) the accuracy of the canonical second order finite difference approx¬ 
imation. 

The observation that the SGF formalism treats the cable exactly then 
leads to the question whether the second order finite difference approximation 
can be derived from it, as both approaches describe the voltage at a given 
location only as a function of the voltage at the neighboring locations and 
the input at that location. We show in a simplified setting that, when the 
distance Ax between the input locations is sufficiently small, the formalism 
reduces to the second order finite difference approximation of the original 
equation. 

Although we expect that the reduction of the SGF formalism to the sec¬ 
ond order finite difference approximation is valid for all equations of the type 
Q and for arbitrary tree graphs, its derivation can become prohibitively 
complex. We therefore constrain ourselves to the passive cable equation ([3]) 
and to the case where this equation is defined on a line of length 1, with a 
homogeneous boundary condition at each end. We suppose that there are 
n — l input locations, distributed evenly with spacing Ax = 1/n. 

\ + StM) ~ v ( x A = Ii(t)S(x - iAx). 

|(0,i) = L#(0,() (50) 

If (l,t) = L B V(l,t) 

Note that the input currents /$(£) can depend on the local voltage V(xi,t), 
as they can be the result of a discretization of the ion channel currents. 
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Second order finite difference approximation To obtain the second 
order finite difference approximation of equation (50), we replace 


d 2 V. V w (t)~2V i (t) + V,. l (t) 

^ - Ax? - 

and we average the inputs for each compartment: 

1 ,(*+1/2) Ax J(t) 

— / umx - iAx)dx = 

Ax 1 / 2 )Ax Ax 


Consequently, we find for % — 1 ,..., n « 1: 

V i+1 (t)-2V i (t) + V i - 1 (t) 


Ax 2 


+ V i (t)-V i (t) = 


m 

Ax ' 


(51) 


(52) 


(53) 


The boundary condition at x = 0 becomes: 

ViO) - V 0 (t) 

- ~Ki -“ La Mt) ’ 

and an analogous expression applies for the boundary condition at x — 1 


(54) 


Reduction of the SGF formalism 

(50) reads: 


The Fourier transform of the problem 


f ® (x,u) - y(u;) 2 F(cu) = E-Li Ii(w)6(x - x { ) 
f (0,ou) = L a (cu)V(0 1 cu) (55) 

[f (1 }U ) = L b (u)V(1, U ), 


where La and Lb are polynomials in u representing the respective Fourier 
transforms of La and Lb , and where 7 2 (oo) = iuj — 1. In the following 
the coordinate u 1 is dropped for notational clarity. We define two linearly 
independent solutions to the homogeneous problem: 


ua(x) = e 7X + kAe 71 
u B (x) = e lx + k B e~ lx 


that satisfy the boundary conditions respectively in x — 0 (ua) and x = 1 
[ub]- Then the Green’s function of problem (55) is given by (Stakgold, 1967 
page 66 ): 


9ij 


UA{xi)u B {x-j) 

-27(fc s -fc^) 

U B (xj)u A (xj) 

-27 {k B -k A ) 


if Xi < Xj 
if Xi > Xj. 


(57) 
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Consequently, the attenuation functions are given by: 


a ij 


U A (xj) 
_ ) U A (Xj) 
U B (Xj) 
U B (Xj) 


if Xi < Xj 
if Xi > Xj. 


( 58 ) 


According to theorem [lj problem (55) is then fully defined by the following 
kernels: 


for i — 0,..., n — 2 

7 _ ^ 2 + 1,2 ^ 2 + 1 , 2+2 ^ 2 + 2,2 

^2+1,2 Z 

-L ^2,2+2 ^2+2,2 

for i = 2,... ,7i — 1 


— 1, 


1 — a i+l,i+2®i+2,i+l 


7 _ ^ 2 , 2+1 ^ 2 , 2 — 1 ^ 2 — 1 , 2+1 

^2,2+1 — \ 7 I ^2,2+1 



1 ^2+1,2—1^2—1,2+1 


1 ^ 2 + 1 , 2 — 1^2 — 1 , 2+1 

hoi = Ooi 

hn,n—1 ®n,n—1 

for i = 1,... ,n — 1 

(1 l®i—l,i) (1 ®i+l,i®i,i+l) 


fi 9r. 


1 l®i—l,i+l 


Using equations (58) and (56), it can be checked that: 


2(k B -k A )smh('y(x j -x i )) jf <- 

■1 _ ) u B (xi)u A (xj) 11 

a 'j°i' ~ ^ ■2{k B -k A )shih^{x l -x 1 )) jf x . > x . 


U A (Xi)u B {Xj) 


Consequently, it follows that: 

• for i — 0,..., n — 2 

• for i — 2,... ,7i — 1 




sinh+Aic) 




sinh( 27 Ax) 

sinh+Ax) 

sinh( 27 Ax) 


(59) 


(60) 


(61) 


(62) 


(63) 


(64) 
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• for i — 1,..., n — 1 


fi = ~ 


1 (sinh^Ax)) 2 
7 sinh( 27 Ax) 


(65) 


With these equations, equation (|12|) for i = 1,..., n — 1 becomes: 

sinh^Ax) 


_ 1 (sinh^Ax)) 2 r | 

vi . ^ * \ J-i \ 


(Vi-i + Vi+i) ■ 


7 sinh(27Ax) sinh(27Ax) 

Since Ax is small, the sinh-functions can be expanded. This gives: 

Ax ~ 1 + ^7 2 Ax 2 


( 66 ) 


Vi 2 + §7 2 Ax 2 Ii + 2+ | 7 2 Ax 2 


(Vi-i + Vi+i) ■ 


(67) 


Multiplying both sides by the denominator 2 + |y 2 Ax 2 , rearranging terms 
and dividing by Ax 2 then gives: 

Vi-i ~ 2 Vi + Vi+i 8 2 ~ 1 2 (fr y \ If 

- 77 Vi + -7 [Vi-! + V i+1 J = —(68) 


Ax 2 


6 


6 


Ax 


Averaging VAi + Vi +1 ~ 2Vi and transforming the resulting equation back to 
the time domain then results precisely in the finite difference approximation 
for % — 1,..., n — 1. 

Let us now investigate equation (12) when i — 0. Here, it holds that: 


Vq — aoiVi. 


(69) 


Substituting the explicit form of a 0 i = m j 4(0)/ua(Ax) in this equation and 
expanding the exponentials up to first order leads to: 


W = 


1 + 


irfc^Ax 


Vi 


(70) 


Since ua satisfies the boundary condition at x = 0, it can be checked that 
7 (1 — fcf0/(l + kjC) = La , so that equation (70) can be rewritten as: 


L a V 0 


Vi-Vq 

Ax 


(71) 


Transforming this equation back to the time domain yields precisely the finite 
difference approximation for Vq. Analogous considerations apply for V n . 
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Implementation 

In this section we summarize the steps one has to undertake to implement 
the SGF formalism. Such an implementation consists of two main parts: an 
initialization part and a simulation part. 

The initialization part consists of the following steps: 


1. For any given set of input locations and a tree graph, a routine has to 
be implemented that can identify the sets of nearest neighbors. Such 
a routine could for instance return a list, where each element is a list 
in itself, containing the input locations that constitute a set of nearest 
neighbors. 


2. The GF’s gij(u) between every pair of elements in the sets of near¬ 
est neighbors have to be computed. We implemented the algorithm of 
(Koch and Poggio, |1985 ) for this purpose. Note that for a set J\f of 
cardinality |jV| , only |j\7| (|A/”| + l)/2 function g i3 have to be computed, 
since g t] = gji (Koch, 1998). From these kernels, the attenuation func¬ 
tions can be derived easily. 


3. The kernels /) and hij have to be computed. This can be done either 
by evaluating the formulas (22) and (23) explicitly (when all matrices 
of which the determinants have to be computed are small) or by the 
derived formulas: 


fi = 9ll A( A/1U...UA/;)- 1 

i4(JVi U...UjV'pU N)~ l 
hj ~ A(N i u ... u Af p u JV)^ 1 ' 


(72) 


Note that in both cases, attenuation function are needed that are not 
directly computed in the previous step, since attenuation matrices of 
unions of sets are considered. This is not a problem, as these attenu¬ 
ation functions can be easily reconstructed from the transitivity prop¬ 
erty. 


4. The partial fraction decomposition of each kernel has to be computed. 
The VF algorithm (Gustavsen and Semiyen 1999[ ) (appendix § |B| is well 
suited for this purpose. When the parameters of the partial fraction 
decomposition are known, the vector Fq, the matrices F\ and Hq and 
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the tensor Hi can be computed, as described in (43) and (44). Since 
Hq resp. H] are sparse in their indices i and j, they can be stored as 
index-number resp. index-array pairs. 

For the simulation part, two main routines have to be called each time step: 


1. A routine that computes the vectors c (t) and d(i) (48). This routine 
also advances synaptic conductances and ion channel state variables 
and its details are well established in the neuroscientific literature (see 


(Rotter and Diesmann 1999) and (Hines, 1984)) 


2. A routine that first computes the vector k(t) (46) and then solves 


equation (49) 


We implemented a prototype of the above outlined initialization algorithm 
in Python. Tree structures were implemented using the btmorph library 


(Torben-Nielsen, 2014). The simulation algorithm was implemented both in 


pure Python and in C++ with a Cython interface. 


3 Validation 


In myelinated axons, the approximation of grouping active currents at a 
discrete set of input locations holds exactly, as the only spots where active 
currents are present are the nodes of Ranvier. These nodes are separated by 
stretches of myelinated fiber of up to 2 mm in length, which can be modeled 
by equation (J3]) . We reproduce the model of (Moore et ah, 1978), where the 
soma, axon initial segment and nodes of Ranvier are equipped with Hodgkin- 


Huxley (Hodgkin and Huxley, 1952) channels, so that the input current is of 
the form: 


Ii(t) = -g Na m i (t) 3 h i (t)(V(x i ,t)-E Na )+g K n i (t) , (V(x i ,t)-E K )-g L (V(x i ,t)-E L ). 

(73) 


Consequently, Ci(t) and d t {t) in equation (48) are given by: 


Ci(t) = g Na mi(tfhi(t)E Na + g K ni(t) A E K + g L E L 
di(t ) = -gN a mi(t) 3 hi(t) - g K ni(t) A - g L , 


(74) 
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Relative parameter f(ms) 


Figure 2: Validation of the SGF formalism on axon (A-B) and den¬ 
drite models (C-F). We implemented the axon model described in (Moore 


et ah, 1978) in the SGF formalism. A: validation of the action potential 


(AP) shape by comparison with an equivalent NEURON model at the axon 
initial segment (AIS) and the 10’th node of Ranvier (NoR). B: reproduction 
of Fig 2 in (Moore et ah, 1978), where the dependence of AP velocity on 


different biophysical parameters is studied. C: The dendritic morphology to¬ 
gether with 50 synaptic input locations (the morphology was retrieved from 
the NeuroMorpho.org repository (Ascoli, 2006) and originally published in 
(Wang et ah, 2002)). D: The number of kernels in the SGF formalism, com¬ 
pared to the number of kernels that would have been required in the normal 
GF formalism. E: The number of operation required per kernel to achieve 
similar levels of accuracy for 3 approaches: Exp, where all the exponentials 
from the VF algorithm are integrated, Quad, where the quadrature is com¬ 
puted explicitly and Mix, where we compute the quadrature for the first K 
(here 3) steps, and use the ODEs to compute the rest of the convolution. 
F: Voltage trace at the soma upon stimulation of the synapses with 5 Hz 
Poisson spike trains. 
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and the variables rrii(t),hi(t) and n^t) evolve according to the following 
equations: 


ihi{t) = 
hi(t ) = 
fii(t) = 


m in f(V(xj,t)) - rrii(t) 
T m (V(Xi,t)) 
hmf(V(Xj,t)) - hj(t) 
T h (V(Xi,t )) 
n in{ (V(xi,t)) — rii(t) 


(75) 


T n (y(xj,f)) 

where m in f, /t in f, n in f, r m , r/,, and r n are functions of the local voltage. 

We implemented this model using the SGF formalism. In Fig [2]A., we 
validate the action potential (AP) shape at the axon initial segment and 
at the 10’th node of Ranvier by comparison with the NEURON simulator 
(Carnevale and Hines, 2006) - the gold standard in neuronal modeling. In 


Fig [2j3 we reproduce Fig 2 in (Moore et ah, 1978), where the dependence of 


the AP velocity on various axonal parameters is investigated. 

Finally, to show that the SGF formalism can also handle complex trees, we 
equipped a stellate cell morphology, comprising an active soma and passive 
dendrites (Fig[2p) (modeled by equation (j3j)) , with 50 conductance based 
synapses, so that 


h{t) = g(Ai{t) ~ Bi(t))(V(xi,t ) - E r ) 


Ait) 




-Aid) + CA'LjSjt - s[ j) ) 
ta 

-b^) + c b e j^d - 4 J) ) 

Tb 


(76) 


where represents the j’th spike time of the input spike train arriving at 
the i’th synapse. In this equation, the variables A and B are used to generate 
a double exponential profile (see for instance (Rotter and Diesmann, 1999[ )). 
In this model, it holds that: 


Cid) = -g(Aid) - Bid))E r 

ddt) = g(Adt) - Bid))- 

Excellent agreement with the NEURON simulator is achieved (Fig [2^). The 
amount of kernels required is far lower than in the normal GF formalism 
(Fig [2j7)) and optimal performance (for which we found ss 7) was achieved 
for K = 3 (Fig@E). 
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Figure 3: The execution time of the SGF formalism compared to 
a neuron simulation. A: The setup, showing the morphology and the 
maximal set of input locations N = 74 input location. Note that this set was 
chosen so that \J\f\ = 2 for all J\f. B: The somatic voltage trace for the first 
second of a 10 s simulation, with a conductance based excitatory synapse at 
each input location indicated in A. C: Execution times as a function of the 
number of input locations. D: The average number of operations per kernel 
as a function of the number of input locations. 


Computational cost 


It is immediately clear that the computational cost of the SGF formalism 


is far lower than the computational cost of the GF formalism (Wybo et al 


2013) from the number of required kernels alone (Fig[2]D). However, the com¬ 
parison of computational cost with the canonical finite difference approaches 
requires a more careful discussion. For n input locations and n t time steps, 
an explicit solver would typically require 0(nkn t n ) steps, whereas a canon¬ 
ical finite difference approach would require 0(d t n t n x ) steps, with n x the 
number of spatial locations at which the finite differences have to be evalu¬ 
ated and dt the maximal degree of temporal differentiation in the operators 
Lf\x). Thus, we expect the SGF formalism to improve performance when 
cnk‘n < dtn x , where c is an a-priori unknown constant that may depend on 
the specific implementation. When the system is stiff, an implicit solver is 
typically needed to guarantee stability. Such solvers require a matrix of size 
n (resp. n x in the case of finite differences) to be inverted each time step, 
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normally an operation of complexity 0(n ^). However, for second order finite 
difference solvers, the tree graph introduces a special structure in this matrix 
(Hines, 1984), so that it can be inverted in 0(n x ) steps. In the SGF formal¬ 
ism, this structure can be maintained if |A/"| = 2 for all J\f and in this case 
the same performance criterion cn^n < d t n x applies as for explicit solvers. 
Note that for higher order finite difference approaches the matrix inversion 
still requires O(n^). 

To test whether this constant c is not excessively large, so that it may 
impede the usefulness of the SGF formalism, we compared the execution 
time of the SGF formalism with the execution time of the NEURON simulator 
(Carnevale and Hines, 2006). In order to obtain a ‘fair’ comparison, attention 
must be given to the exact type of implementation. For instance, comparing 
an explicit simulation paradigm with NEURON makes little sense, as NEURON 
implements a semi-implicit paradigm and would be severely at the disadvan¬ 
tage. On the other hand, NEURON makes use of the Hines algorithm to invert 
the systems’ matrix (Hines, 1984), and distributing input location completely 
at random in the SGF formalism (and using a semi-implicit paradigm as well) 
would exclude the use of this algorithm, which would severely disadvantage 
the SGF formalism. We therefor opted to implement a simulator in C++, 
that implements the same semi-implicit method described in equations (49) 
and (47), but where the inversion uses the Hines algorithm. This restricts the 
input locations to configurations where \J\f\ = 2 for all A f. We distributed 2 
to 74 input locations in such a configuration (Fig[3|\ shows this configuration 
for n = 74) on the same stellate cell morphology used in Fig [2] (panels C-E). 
No active channels where added to the model. In NEURON, the cell model 
contained 299 compartments for a total dendritic length of approximately 
4000 /irn ( Ax ~ 13.5 /iin). We added one excitatory, conductance based 
synapse at each input location and gave it a Poisson spike train of a certain 
rate, so that the total presynaptic rate at all synapses together was 1000 Hz. 
In Fig [3^3 we show the first second of a 10 s simulation where n = 74. It can 
be seen that both traces agree impeccably. We then compared the execution 
times of both models for an increasing number of input locations (Fig[3p) 
for a simulation time of 10 s at a time step of 0.1 ms. It can be seen that 
the SGF implementation outperforms the NEURON model until n ~ 55 (on a 
machine equipped with an Intel Core i7-3770 CPU @ 3.4GHz and 16 GB of 
RAM, running ubuntu 14.04 LTS). In Fig[3p, we plotted (n^), the average 
number of operations per kernel and per time step. It can be seen that this 
number correlates with the execution time. When there are few input loca- 













tions, kernels are in general longer and hence this number is higher. Then, 
when more kernels are added, this number drops quite rapidly. This explains 
the small slope of the execution times in the SGF formalism until n ~ 35, 
associated with the rapid decrease of (n^). After that, (n^) decreases slower, 
and hence the slope of the execution times becomes steeper. 


4 Summarizing remarks 


We have proven that on tree graphs, the GF formalism for linear, non- 
homogeneous, time-invariant PDE’s with inputs at a discrete set of n well- 
chosen locations requires only 0(n) rather than n 2 kernels, and termed this 
the SGF formalism. We discussed the meaning of ‘well-chosen’, namely that 
the sizes of the sets of nearest neighbors must be small (ideally |A/"| = 2 for 
all A/"). We have shown furthermore for equation (J3]) that in the limit of 
a small distances between the input locations, the SGF formalism can be 
reduced to the second order finite difference approximation. Thus, in some 
sense, the SGF formalism can be seen as a generalization of the second order 
finite difference approximation to arbitrary distance steps (as long as what 
lies between the input locations is approximately linear). We also employed 


the VF algorithm (Gustavsen and Semiyen 1999), that fits the kernels with 


sums of exponentials, to design an efficient simulation algorithm. 

We then validated our algorithm on two neuroscientific problems: the 
modeling of myelinated axons and of dendritic integration. We showed that 
there was excellent agreement between the SGF formalism and the de facto 
standard NEURON simulator. We found furthermore that, even for complex, 
branched morphologies, excellent accuracy was achieved for < 15. We 
also discussed when an increase in computational performance is expected 
from the SGF formalism, and furthered our findings by comparing an efficient 
C++ implementation of the SGF formalism with the same model in NEURON. 
Finally, we discuss in a broader sense when the SGF formalism is expected 
to provide advantage over other approaches. 


5 Discussion 

It can be seen that the question whether one should prefer the SGF formalism 
over second order finite difference approaches depends on the values of Uk 
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and n, which are not necessarily independent. Fig [3] illustrates this: for the 
same dendritic arborization, and using the passive cable equation (dt = 1), 
rik ranges from « 14 for n = 2 to ~ 4 for n = 74. This leads to n^n ps 28 to 
296 operations per time step. The dendritic arborization on the other hand 
has a total length of approximately 4000 /mi. According to the Lambda rule 
(Carnevale and Hines, 2006), the canonically accepted method of calculating 
the distance step between compartments, Ax ~ 10 — 15 /mi. Consequently, a 
second order finite difference approximation would use approximately n x ps 
300. Hence this estimate would indicate that performance can be gained 
when n < 75. That we found this number to be slightly lower (n < 55) 
in the simulations we conducted (Fig |3j) may be due to implementation and 
optimization related factors. One important remark is that for efficient semi- 
implicit simulations, the input locations must be chosen so that \J\f\ = 2 
for all J\f to be able to use the Hines algorithm. The question may now 
be asked whether this is overly restrictive. We believe that this is not the 
case: placing an input location at the soma for instance separates all main 
dendrites. Bifurcations in the dendrites may induce further sets of nearest 
neighbors of sizes slightly larger than 2. Whenever this occurs however, an 
input location can be added at the bifurcation point, so that such a set is 
split as well. 

Thus the choice of the SGF formalism over the finite difference approx¬ 
imation depends on the neural system at hand. In some systems, such as 
bipolar neurons used in auditory coincidence detection (Agmon-Snir et ah] 


1998 Wybo et ah, 2013), inputs occur only at a small set of locations, and 
the cells’ computation is performed by linear membrane properties. Other 
systems, such as myelinated axons, have non-linear membrane currents that 
are concentrated at a discrete set of locations - the nodes of Ranvier, with 
stretches of myelinated fiber in between that behave approximately linear. 
It is to be expected that in such systems, the SGF formalism may signifi¬ 
cantly improve performance over the second order finite difference approach. 
In other systems, such as cortical cells, inputs are distributed throughout 
the dendritic arborization in an almost continuous fashion. In these cells, 
the answer to the question whether the SGF formalism yields computational 
advantage is negative, when one aims at retaining all biophysical detail. In¬ 
deed, it is clear that when the number of input locations in the SGF formal¬ 
ism approaches the number of compartments NEURON requires, there is no 
computational gain in using the former. Nevertheless, computational neu¬ 
roscientists have been seeking ways of reducing the cost of simulating these 
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cells to be able to use them in large scale network simulations. Most often, 
they aim to achieve this by drastically reducing the number of compartments 
(Traub et ah, 2005). This approach however also changes the spatio-temporal 
response properties of the nerve cells. With the SGF formalism, inputs that 
would otherwise be grouped in a small number of compartments may now 
be grouped at a small number of input locations, while the response proper¬ 
ties induced by the neuronal morphology would remain unchanged. Finally, 
in dendritic arborizations, a number of resonance phenomena have been ob¬ 
served that can be modeled with linearized (quasi-active) ion channels ( Koch[ 
1984 Laudanski et al.[ 2014). Since the transitivity property (see appendix 
§|Bj)) holds for equations of this type as well (of the general form (J8])) , such 
systems can be modeled implicitly in the SGF formalism. Nevertheless, the 
validity of this linearization depends on the size of the fluctuations around 
the operating point and the ion channel under study and has to be checked. 

In the SGF formalism there is a significant initialization phase before a 
neuron model can be simulated (see Implementation). The computational 
cost of this phase is higher than the cost of the initialization phase for finite 
difference approaches, and depends on the complexity of the tree graph and 
the number of input locations. We thus expect that our SGF formalism will 
be advantageous in use cases where frequent re-initialization of the model 
is not required. This is typically the case for network simulations, where 
a limited number of prototypical nerve cells may be initialized and used 
throughout the network. 

Another important matter, next to computational cost, is the accuracy 
of the SGF formalism. While the sparsihcation is exact, the transform back 
to the time-domain along with the specific integration algorithm might in¬ 
troduce errors. The use of the approximate VF algorithm however impedes 
a systematic analysis of these errors. Nevertheless, we found that typically 
this error was very low (Fig [4^3, C). Furthermore, after application of the 
integration paradigm described in this work, where the convolutions with 
exponentials are integrated analytically assuming that the voltage varies lin¬ 
early in between grid points, we found that all our numerical experiments 
agreed very well with equivalent NEURON simulations. 


31 












Acknowledgments 

We thank Luc Guyot for many insightful comments on the manuscript. The 
research leading to these results has received funding from the European 
Union - Seventh Framework Programme [FP7/2007-2013] under grant agree¬ 
ments no. 269921 (BrainScaleS) and no. 604102 (The HUMAN BRAIN 
PROJECT) as well as funding from the EPFL Blue Brain Project Fund and 
the ETH Board to the Blue Brain Project. 


A Proof of Lemma [T] 


Proof. The PDE defined on a line 

Consider a line of length L (0 < x < L), which trivially has two leafs 
(A = {Ai, A 2 }). Fourier transforming a PDE such as (J8| leads to a boundary 
value problem of the form: 


~ <9 2 ~. ~ d ~ ~ 

L 0 (x,lo)^—V(x,co)+L 1 (x,co)—V(x,uj)+L 2 (x,co)V(x,uj) = I{x,u), 0 < x 

ox z ox 

(78) 

B x 'V(0,u) ■= Ll‘(ui)f-V(L,w) + Li l (u)V(0,u) = I x ‘(ui) 

% x ' (79) 

B Xx V( 0,w) := V(L,u) + L^(u)V(L,u) = /*», 

The Green’s functions g(x,Xo,oj) is obtained from solving this problem for 


I(x, u) = 5(x — x 0 ) and I Xi (lo) = /A 2 (u) = 0, and is given in ( |StakgoldHl967| 
page 66) : 

f U 1 {x)u2(x 0 ) D < T < Tr, 

n ( T T \ — ) a 0 (x 0 )W[u 1 ,U2\x 0 ) ’ — — 0 ('on) 

9{x, xq) — uUxn)M Xq < x < Li 


ao( k xo)W(ui,ur,xo )' 


where Ui(x) is a non-trivial solution of the homogeneous equation satisfying 
B Xl u\ = 0 and w 2 (a;) a non-trivial solution satisfying R A2 u 2 = 0, where 
W(ui,u 2 ] x) denotes the Wronskian of both solutions evaluated at x and 
where we have omitted the dependence on cn for clarity. From this equation, 
it can be checked that the transitivity property holds: 


g(x!,x 3 ) 


g(x 1 ,x 2 )g(x 2 ,x 3 ) 
g(x 2 ,x 2 ) 


(81) 


< L 
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when X]_ < x 2 < x 3 or x% < x 2 < X\. 


The generalization to a tree graph 


Consider now the generalization of problem (78) to a tree graph. On each 
edge e G E, an operator of the form: 

~ , d 2 


dx 2 


constrains the field V e : 


where it is understood that x signifies the space coordinate on the edge under 



(82) 

I e (x), 

(83) 


consideration. On each leaf A G A a boundary condition of the form (79) 
holds: 

B X V X = I x (84) 

and at each node that is not a leaf (ft G <f>: 


V e = V e ', Ve, e' G E{(ft) 

£«*<*> = i*, 


(85) 

( 86 ) 


Equation (85) assures continuity of the field and a condition of the form (86) 


is imposed in many physical systems to assure the conservation of flow. 

We wish to determine Green’s function g(x,x o) of this problem. When 
xo is located on edge eo, we need to solve this problem for I e (x) = 5 eeo 5(x — 
xq ),/ A = 7^ = 0. We will first show that at each end of edge eo boundary 


conditions of the form (79) hold, by using (84), (85) and (86). 


To do so, we only need the following recursion rule: Consider a node 
(ft, and suppose that all but one of the edges in E((ft) satisfy a boundary 


condition B e V e = 0 of the type (79) at the opposite end of node (ft. Within 
each edge, we chose the spatial coordinate x e to be 7/ (i.e. the edge’s length) 
at that far end and 0 at the node. Let us call the edge that does not satisfy 
the boundary condition e'. Thus we have: 

B e V e = L\^-V e (L e ) + ^W(77) = 0, Ve G E{(ft) \ e', (87) 


and from (85) and (86) if follows that: 


W'(0) = W(0), Ve G E((ft) 

0) + £ e '£v'(0)=0, 
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( 88 ) 

(89) 














We will show that from conditions ( 88 ) and (89) a boundary condition of the 
form (79) can be derived for V £ (0) (i.e. the held on edge e' at the location of 


node 0), as long as the differential equations on edges e G E((j))\e' are homo¬ 


geneous. Let u e (x) be a non-trivial solution of the homogeneous problem (83) 
on edge e that satisfies condition ( |87| ). Every solution on that edge is neces¬ 
sarily of the form V e (x ) = A e u e (x). As a consequence of ( 88 ), A e = -—— 


which leads to a constraint on the derivative yW e ( 0 ) = 


equation (89) becomes: 


v e (o) a 

« e (0) dx 



V‘‘(0) + l-Av‘‘(0) = 0, 


M e (0) ’ 

M e (0). Thus, 


(90) 


precisely the boundary condition for V e> we were after. 

Applying this operation recursively throughout the tree, starting from 
the leafs until edge e 0 , assures that this edge has a boundary condition of the 


form (79) at both ends. Thus, on this edge, g(x,x 0 ) is of the form (80). 


To prove the transitivity property (81) for two arbitrary points X\ and £ 3 


on the tree graph, and for a point X 2 that is on the shortest path between x\ 
and £ 3 , we distinguish four cases. 

1. £1 , £ 2 , £3 are on the same edge. Since the segment has boundary 


conditions of the type (79), the Green’s function may be constructed 


as in (80), and thus (81) holds 


£1 and £2 on the same edge, £3 is on a different edge. Let 0 be 

the node adjacent to the edge eo on which £1 and £2 are located and on 
the shortest path between £2 and £ 3 . Necessarily, the Green’s function 
at that point satisfies 


0(0, aq) = 


g(^,x 2 )g(x2,x 1 ) 

g(x 2,x 2 ) 


(91) 


which then determines the solution on the adjacent edges e G E((ft)\e 0 
(where we choose the spatial coordinate x e to be 0 at 0 and L e at the 
opposite end). On either of these edges, the solution is of the form 
I, with u e (x £ ) a solution satisfying the derived condition of type 


A £ u e (x e 


(79) at the opposite end. Condition (85) then imposes A £ = = 
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g(<M2)gQE2,a:i) J 
g(x 2 ,x 2 )u*(0) ’ allCl tllUS ' 


/ e x 0(0, ah) 0(0,®2)0(®2,a;i) e , £ x 

g[X ,Xi) = - ——'u (x ) = — -r ————u (x ) 


u e ( 0 ) 

•,X 2 )u‘ 
« e (0) 


g(x 2 ,x 2 )u e (0) 


g(0,.i 2 )u (x lg( X2 ^ Xl ) g(x e ,X 2 )g(x 2 ,Xi) 


(92) 


g(x 2 ,x 2 ) g(x 2 ,x 2 ) 

We may apply this consideration recursively through the tree graph, 


until we arrive at the point X 3 , which proves relation (81) in this case. 


3. x 2 and £3 on the same edge, x\ is on a different edge. Let e 

denote the edge on which x 2 and £3 are located, and let us denote 
the node adjacent to that edge at the side of x\ by 0 , and take the x- 
coordinate describing the position in that edge to be zero there ( x e = 0 ). 
Then W(0) = g(0, X\). On the other side of the edge, at x e = L € , the 


derived condition of the form (|79|) holds, and thus the field in that edge 

_ g( 0 ,S 1 ), 


satisfies V e (x e ) = 


u e (0) 


-U 


(£ e )(= g(£ e ,£i)), where u e (x e ) is a solution 


satisfying this boundary condition. The Green’s functions g(x 2 , x 2 ) and 


g(x 3 ,x 2 ) are still of the form (80), as in this case the derived boundary 
conditions are valid on both ends of the edge. Let v e (x e ) be a solution 
that satisfies the homogeneous boundary condition at £ e = 0. Then, 


using (80), the following holds: 


g{x 3 , X 2 ) — - -g(x 2 , £ 1 ) 

g{x 2 ,x 2 ) 


v £ (x 2 )u € (x 3 ) a 0 (x 2 )W(u £ , v e , £ 2 ) g{4>, £i)u e (£ 2 ) 


a 0 (x 2 )W(u e ,v e ,x 2 ) v e (x 2 )u e (x 2 ) 

= 0(0,Xi)m 6 (£ 3 ) 

u e (0) 

= g(x 3 ,xi) 


u e (0) 


(93) 


4. £ 1 , £ 2 , £3 are on different edges. This case is proved by combining 
the two previous cases. 
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B Vector fitting 


Here we briefly explain the VF algorithm (Gustavsen and Semiyen 1999) as 
we implemented it. The version of this algorithm we needed approximates a 
complex function f(s), for which |/(s)| —> 0 when | .s\ —> oo and|/(s)| > 0 , as 
follows: 

oil + s’ 


/(«> 


(94) 


where the parameter L is chosen. It does so in two steps: First the poles cp 
are identified and then the residues 7 / are fitted. 


Pole identification 

First, a set of chosen starting poles cp is specified, and an unknown auxiliary 
function er(s) with these poles is proposed, so that: 


ff ( s ) = E 

1=1 

H s )f( s )\ At = 


11 


Oil 
L 

E- 


+ 1 


ii 


l = l Oil + s 


(95) 

(96) 


Multiplying equation (95) with /(s), and equating this with equation (96) 
gives: 


V- 1 -/ v- f( s ) - t! ^ 

2 ^ 1 ~ 2^ =——■H = f(s). 

l = l Oil + s J =1 Oil + s 


(97) 


When / is evaluated at enough frequency points s (we found it sufficient 
to sample /(s) on the imaginary axis s = 6 l,j = 1,..., N on an 

equidistant scale for small u 1 and on a logarithmic scale for large 00 ), this 
gives an over-determined system that can be solved for 7 1 and 7 ^ by the 
least squares method. The poles of f(s) are then given by the zeros of a(s) 


since, as the parametrization for cr(s) satisfies equation (£ 

7), it holds that 

(see (Gustavsen and Semiyen 

1999 Hendrickx and Dliaene, 

2006 

) for further 


details): 


f(s) 


W(s)f{s)] fit 

a(s) 


(98) 
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Figure 4: Illustration of the VF algorithm. A: A typical kernel is accu¬ 
rately approximated by the VF algorithm (here L = 10). B: The error E(oj) 
of this fit. C: The average error of the fit as a function of L. 


These zeros can be found as the eigenvalues of the matrix: 


H = 


( 


V 


«1 - 71 “72 

—71 02 — 72 


~1L 
-7 L 


\ 


-7l 


-72 


OIL 


(99) 


7 L 


Note that this procedure can be part of an iterative optimization, where the 
newly found poles can be used as the starting poles for the next iteration. 


Residue fitting 

Once the poles cq are known, the residues can be determined by solving the 
over-determined system: 


^—7 1 = f(s) (100) 

+ s 

for 7 1 by the least squares method. 

Accuracy 

The VF algorithm does not provide the number of poles for the fit, nor does 
it provide an estimate for the accuracy of the fit with a given number of 
poles. We chose the number of poles as the smallest number for which the 
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approximation gave a sufficient accuracy, defined as: 


E(ujj) 



- < e, j — 1,..., N. 

maXj f(iujj ) 


(101) 


For our nerve models, we found that e = 10“ 8 was sufficient. Usually, this 
accuracy was reached with 10 < L < 20. In Fig [I]A, we show a typical 
example of a kernel from the SGF formalism, together with its approximation 
with L — 10. In Fig[4]B we show the error of this approximation as a function 
of the frequency, whereas Fig[4p we show how the average error, defined as: 



( 102 ) 


changes with L. 
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